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Abstract 


Accurate, reliable and efficient simulations of moisture fluxes through porous 
media are desirable in hydrological and environmental studies, as v/ell as in civil and 
environmental engineering. For these simulations we need reliable estimate of flow 
parameters. Several mathematical models are available, which attempt to describe various 
soil-hydraulic parameters. Among these parameters, hydraulic conductivity K is the one, 
which still needs a serious attention. The most acceptable mathematical prediction of K is 
given by a model, originally proposed by Mualem and then modified by van Genuchten 
and hence known as VGM model. The correct prediction of K depends on the model 
parameters a and n and hence the suitable values of these parameters are required to be 
determined first. This dissertation is aimed at determination of suitable values of a and n 
from column outflow experiment. The outflow is simulated for a possible range of a and 
n and then by the analysis of these outflow patterns two quantities, namely, T50 and slope 
($60-40), have been foxmd to have a coherent relationship with VGM parameters, hence, 
are suitable for calibration of these parameters. A graphical representation has been given 
from which value of a and n can be directly interpreted. 



Chapter 1 


Introduction 


1.1 General 

The desire and necessity to explore nature for our civilizational development 
ignites the zeal to understand the enigmatic style, in which water, one of the basic 
element of nature, behaves. Moreover, if the water is in form of groundwater, the 
phenomenon becomes highly erratic. 

With increasing demands on groundwater resources the need for an accurate 
prediction of the subsurface flow and chemical species transport under different hydro- 
geological, climatic and ecological conditions have greatly accentuated the need to 
understand these processes and to evaluate effects of management practices and 
remediation techniques. Accurate, reliable and efficient simulations of moisture fluxes 
through porous media are desirable in hydrological and environmental studies, as well as 
in civil and environmental engineering. The ability to model time dependent flows in 
composite soil formations that may be intermittently saturated and drained is particularly 
important from the point of view of physical realism. Several computer codes based on 
numerical models have been used, but with these powerful tools, comes the need for the 
ability of accurate determination of required model parameters. In some cases, only 
groundwater flow is of interest, and the saturated hydraulic conductivity, Ks, must be 
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determined. However, the processes of flow and transport in the vadose zone are essential 
because most groundwater contaminant sources originate in this zone. Model calibration 
in the unsaturated zone may be particularly difficult due to problems in formulating the 
constitutive relations for this special type of two-phase flow, namely the water retention 
curve and the unsaturated hydraulic conductivity. 


1.2 Techniques to determine relationship between h, ^and K 

There are some laboratory and field methods to determine the relationships 
between the pressure head h, the water content 0, and the hydraulic conductivity K. 
Traditionally, direct steady-state methods for the determination of these highly non-linear 
functions exist [Klute and Dirksen (1986)], but recently, transient experimental methods 
coupled with inverse modelling techniques have become more attractive [Kool et al 
(1987), Yeh (1986)]. This parameter identification technique involves the numerical 
solution of the water flow equation for unsaturated/saturated porous media, subject to the 
imposed initial and boundary conditions. The constitutive relations, the so-called 
hydraulic properties are assumed to be described by analytical functions characterized by 
a limited number of parameters. During an experiment some auxiliary variables are 
measured, e.g. cumulative outflow, pressure head, water content, or infiltration. Then, the 
priori unknown parameters are determined by minimizing the objective function 
containing the deviations between observed and predicted quantities. 

The use of laboratory outflow experiments for estimation of unsaturated hydraulic 
properties is advantageous because it is flexible in initial and boundary conditions and 


2 



not very time-consuming. One-step outflow experiments (OS), where an initial saturated 
soil column is drained by a one-step pressure change at the lower boundary, were used in 
conjunction with inverse modelling techniques first by Parker et al. (1985) and Kool et al 
(1985). If only cumulative water outflow is measured there, the inverse problem could be 
ill-posed and lead to a non-unique solution {Zurmuhl (1996)]. Finally, it appears that 
either equilibrium data are needed {Van Dam et al (1992)], or measurements of pressure 
head at one point inside the soil column can improve the performance of the OS-method 
[Toorman et al (1992)]. However, OS methods on large columns can produce dynamic 
capillary pressure-saturation relationships depending on the lower boundary value 
changes [Vachaud et al (1972), Staffer (1978), Nutzmann et al (1994)]. The reasons for 
that and a theoretical model of such relationships were discussed by Hassanizadeh 
(1997). 

These considerations stimulated the investigation of the multi-step outflow 
method (MS), which uses small pressure steps to induce drainage of the soil column. 
Superiority of the MS to the OS-method on small columns was reported in Van Dam et al 
(1994). They showed that the MS experiments with only cumulative outflow data contain 
sufficient information for the unique determination of the soil hydraulic functions, using 
initial estimates derived from the outflow experiment itself. Eching and Hopmans (1993) 
found that the inverse solution technique is greatly improved when cumulative outflow 
data are supplemented with simultaneously measured pressure head data from some 
position inside the column during the MS-expeiiment 

To study the reason for this behaviour the mathematics of inverse problems has to 
be considered. Carrera and Neuman (1986) defined criteria of identifiability and 
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uniqueness according to which inverse problems are well-posed for determining the 
aquifer parameters of groundwater flow. Russo et al (1991) and Toorman et al (1992) 
analyzed conditions of well-posedness by evaluating the response surfaces of an inverse 
problem to estimate the unsaturated hydraulic properties. As shown by Mous (1993), the 
non-uniqueness of the estimates is not due to a bad choice of the optimization algorithm, 
but is merely a consequence of the structure of the model and the design of the 
e.xperiment. Based on a rank-analysis of the Hessian matrix conditions of local 
identifiability could be proved, and the number of identifiable parameters related to the 
experiment was calculated. Zurmuhl (1996) investigated parameter identifiability and 
uniqueness on OS and MS experiments with respect to sensitivity coefficients, and 
showed that only the MS-method can produce uncorrelated and thus linear independent 
parameters. 

One-step and multi-step experiments commonly were carried out on small in-situ 
soil samples. Eching and Hopmans (1993) pointed out that the optimized soil hydraulic 
functions as determined firom soil cores do not necessarily represent in-situ soil 
behaviotir. This may be due to the heterogeneity of undisturbed soil cores and their small 
sizes [Kool sad Parker (1988)]. 


1.3 Mathematical models and Richard’s equation 

Mathematical models typically use the Richard’s equation to describe variably 
saturated flows. It is defined by coupling a statement of flow continuity with the Darcy’s 
equation and is commonly cast into one of the following forms: 
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dK{e) 

dz 


(Moisture-based 6i-form) 


(I.l) 
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dz 


Die) 


dz 


c(;i)*=l. 

dt dz 


Kih) 


m 

dz 


dKjh) 

dz 


dt 


dz 


KiO) 


dz 


dKjh) 

dz 


(Pressure-based, h- form) (1.2) 


(Mixed or coupled form) (1.3) 


where h is the pressure head [L], 9 {h) is the volumetric moisture content, t is time 
[T], z is the (positive downward) depth [L], Kih) or Kiff) [L/T] is the unsaturated 
hydraulic conductivity, Cih) = d^d^ [1/L] is the moisture capacity and Di9)=Ki6)/Ci9) 
is unsaturated diffusivity [L^/T]. 

The solution of Richard’s equation requires the specification of soil constitutive 
functions: 

• Hydraulic conductivity KQi) [L/T], 

• Specific capacity Cih) = dOIdh [L''] and 

• Diffusivity = Ki9)ICie) [L^/T]. 

The combination of highly non-linear constitutive functions with non-trivial 
boundary and initial conditions precludes all but the most simplified analytic approaches 
to the solution of (1.1), (1.2) and (1.3). Common practical approaches for the analysis of 
variably saturated flows are mixed and pressure-based numerical formulations, which 
employ low-order finite difference or finite element spatial discretisation and simple 
Euler time stepping {Celia et al (1990), Paniconi et al (1991), Rathfelder and Abriola 
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(1994), Bergamaschi and Putti et al (1999)]. The numerical stability of the temporal 
approximation is enforced by employing implicit (typically, backward Euler) time 
stepping, while oscillations in the finite element spatial discretisation are controlled by 
lumping [Celia et al (1990), Ju and Kung (1997)]. Substantial research has also been 
dedicated to the solution of the non-linear discrete systems that arise in implicit time 
stepping schemes [Paniconi et. al. (1991), Lehmann and Ackerer (1998), Miller et al 
(1998)], and to the incorporation of the soil constitutive relationships [Miller et al 
(1998)]. 

An appreciable advantage of a numerical scheme based on the mixed form of 
Richard’s equation is its inherent conservation of mass. Conversely, standard numerical 
approximations based on the pressure form of Richard’s equation develop undesirable 
mass balance problems [Celia et. al. (1990)], seriously undermining their physical basis. 

In practical applications where the computational speed of the solution is an 
important priority, numerical errors of the order of 0.1-1% are typically acceptable. In 
these cases low-order schemes are competitive with higher-order solvers, and may in fact 
be preferred due to their better stability and algorithmic simplicity [Wood (1990)]. 
Numerical experiments have shown, for example, that the second-order accurate 
Crank-Nicolson scheme outperformed the first-order accurate backward Euler scheme 
only when relative errors below 0.005% were required [Wood (1990)]. The results of 
other numerical investigations [e.g., Paniconi et. al. (1991) and Toed et. al. (1997)] also 
suggest that low-order schemes are competitive with higher-order schemes when coarse 
time steps are used. Indeed, most practical codes implement simple first- or second-order 
^proximations [Celia etal (1990)]. 
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Chapter 2 


Outline of present problem and 

LITERATURE REVIEW 


2.1 Introduction 


For 1-D vertical flow, Richard’s equation is given by: 


de{K) 

dt 


_5 

dz 


K(h)^-K(h) 

OZ 


( 2 . 1 ) 


Water movement through the unsaturated zone is commonly analyzed by solving 
Richard’s equation {Richards (193.1:)]; Analytical and simplified solutions- of Richard’s: 
equation [Philip (1969), Parlange (1972), Broadbridge et al (1988), Warrick (1991)] 
provide useful tools for studying simple unsaturated flow systems with relatively simple 
initial and boimdary conditions. The solutions of these models are based on the following 
assumptions: 

(i) Soil is homogeneous, 

(ii) Initial moisture content is uniform throughout the soil profile, and 

(iii) Moisture content at the soil surface is constant and rainfall or irrigation 
rate is constant. 

In addition, models give accurate results only for a particular type of soil. For 
example, Green-Ampt Model, which is based on the assumption of saturated plug flow, 
fails in situations where a coarse textured soil with high hydraulic conductivity underlies 
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a fine textured soil with low hydraulic conductivity [Chow (1988)]. In the field, soils are 
seldom homogeneous, initial moisture content is seldom uniformly distributed and in 
most field situations during rainfall or irrigation, the soil surface is rarely at constant 
saturation. For accurate prediction of moisture movement under realistic boundary 
conditions in field soils, one has to resort to numerical models, which are versatile in 
handling the non-homogeneity and different kinds of boundary conditions. 

Then comes the measurement of the unsaturated conductivity. Accurate 
measurement is generally cumbersome, costly, and very time-consuming. Consequently, 
many attempts have been made to develop indirect methods that predict the conductivity 
function from the more easily measured water retention curve [Mualent (1976), Leij 
(1996)]. Most of the predictive conductivity models are based on the assumption of 
having an ideal capillary medium characterized by a certain pore-size distribution 
function. Although not necessarily valid for all soils, this assumption is widely accepted 
as an effective working hypothesis [M/a/ew (1976)]. 

A variety of empirical equations have been used to describe the soil water 
retention curve. One of the most popular equations is the power-law function initially 
proposed by Brooks and Corey (1964). This equation leads to an air-entry value, ha, in 
the soil water retention curve above which the soil is assumed to be saturated. The 
Brooks and Corey equation was modified by van Genuchten (1980) to enable a more 
accurate description of observed soil hydraulic data near saturation, especially for 
undisturbed and many fime-textured soils. The choice of the analytical model for 6(h') can 
significantly affect the predicted K(h) function obtained with one of the statistical pore- 
size distribution models. One reason for this is that the predicted conductivity function is 
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extremely sensitive to small changes in the shape of the retention curve near saturation 
[Vogel et. al. (2001)]. This sensitivity is a major cause of the sometimes significant 
differences between predicted K{h) functions obtained with the Brooks and Corey and 
van Genuchten retention functions. The differences are specially important for fine- 
textured soils which can exhibit extreme non-linearity in Kill) close to saturation when 
Van Genuchten’s equations are used. The differences are generally much less severe for 
coarse-textured soils. 

The presence of highly non-linear K{h) relationships near saturation can also 
substantially impact the performance of numerical solutions of Eq. (2.1) in terms of the 
accuracy, stability, and rate of convergence of the invoked numerical scheme [Vogel et. 
al. (2001)]. While numerical solutions for solving the variably saturated flow equation 
have been significantly improved in recent years [Milly (1985), Celia et. al. (1990), 
Kirkland et. al. (1992), Huang et. al. (1996)], most improvements focused on numerical 
problems associated with the infiltration of water in very dry coarse-textured soils [hills 
et. al. {19Z9), Huang et. al. (1996)]. 


2.2 Purpose of present study 

Several numerical models have been developed for simulating water movement in 
unsaturated porous media using finite difference, finite element, and integrated finite 
difference methods. Most of the numerical models considered focus their attention either 
on improving the existing methods or on concentrating on one process such as 
infiltration, gravity drainage or evaporation. However, very few attempts have been made 
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to study the sensitivity of different processes with respect to the unsaturated soil 
parameters. The objective of the present study is 

• to analyze VGM model and study the effect of model parameters (a and n) 
on 1-D transient flow in vadose zone. 

• To develop a methodology for identification of the VDM parameters for 
given outflow data. 


2.3 VGM Model 


The predictive model of Mualem (1976) for the relative hydraulic conductivity 
function, Krih), may be written in the form: 


Kr (S,) = Si 


1 

1 

7 

1 

. 0 

/ 



( 2 . 2 ) 


where = K/Ks, K is the unsaturated hydraulic conductivity, Ks the saturated 
hydraulic conductivity, / the pore connectivity parameter usually assumed to be 0.5 
following AfMfl/e/M (1976), and Se is the effective saturation given by 


' e-e. 


(2.3) 


where 9 Qi) is the volumetric water content, and 9 and 9^ are the residual and 
saturated water contents, respectively. Substituting van Genuchten's (1980) expression 
for the soil water retention curve, i.e.. 


9{h). 


9 .^ 


9 - 9 . 


\ + {a-hy 


9, 


h<0 

h>Q 


(2.4) 
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into Eq. (2.2) leads to the following equation for the unsaturated hydraulic 


conductivity 


K{h) = 


K,K^{h), 

K, 


h<0 

h>0 


(2.5) 


where 

^.(5.) = 5'[l-(l-5r'")'"]^ (2.6) 

in which n and a are empirical shape parameters (with « > 1), and /n = 1 - \/n. 
The restriction m = 1 - l/« is necessary to permit direct integration of Eq. (2.2) once the 
inverse of Eq. (2.4) is substituted into Eq. (2.2). Much more complicated expressions for 
K^{h) result when the parameters m and n are assumed to be mutually independent [Vogel 
et al (2001)]. In this dissertation, Eqs. (2.4) to (2.6) will be referred as the VGM model. 


2.3.1 Shape of hydraulic function 

Now let us see how VGM equations describe the soil hydraulic functions. Fig. 
(2.1) to (2.3) show variations in the soil hydraulic characteristics for a series of «- values, 
while keeping all other parameters constant (a = 0.005 cm'*, 0s = 0.40 and 0r = 0. 10). 

Notice that the K^h) flmction in Fig.2.3 exhibits an abrupt drop at saturation when 
n becomes less than about 1.5. The magnitude of the decrease at = 0 depends not only 
on the parameter n, but to a lesser extent also on cc, which may be viewed as a scaling 
factor for h. The extreme non-linearity occurs only when \< n <1. Close inspection of 
Fig. ( 2.2) and (2.3) shows that the soil water capacity function, C{h)=d.6 /6h, and relative 
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Pressure Head, (cm) 


Fig. 2.3 Plot of Kt vs. h 


conductivity function, Kr(/i), both change their shape near saturation when n changes 
from n < 2 to « > 2. Vogel et al (2001) has discussed that the slopes of these two 
functions change from -oc for C(h) and oc for KrQi) when « < 2 to some non-zero finite 
values when « = 2, and to zero when n>2. This effect of the value of n on the shape of 
the retention curve near saturation was discussed previously by Luckner et al (1989). 
They showed that the restriction « > 1 guarantees only first-order continuity in 9 Qi) at 
saturation, while n> 2 ensiures also second-order continuity in 6{h), and hence first-order 
continuity in CQi). 
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2.3.2 Values of parameters or and n 


As the soil properties tend more towards sand, the value of a and n increases. 
Average values of VGM soil hydraulic parameters for 12 major soil textural groups, 
according to Carsel and Parrish (1988) are given below: 


Table 2.1: 


Texture 

e, (mVm^) 

(9s (mVm^) 

a (cm’') 

n 

ATs (cm/day) 

(1) Sand 

0.045 

0.430 

0.145 

2.68 

712.80 

(2) Loamy sand 

0.057 

0.410 

0.124 

2.28 

350.20 

(3) Sandy Loam 

0.065 

0.410 

0.075 

1.89 

106.10 

(4) Loam 

0.078 

0.430 

0.036 

1.56 

24.96 

(5) Silt 

0.034 

0.460 

0.016 

1.37 

6.00 

(6) Silt loam 

0.067 

0.450 

0.020 

1.41 

10.8 

(7) Sandy clay loam 0. 1 0 

0.390 

0.059 

1.48 

31.44 

(8) Clay loam 

0.095 

0.410 

0.019 

1.31 

6.14 

(9) Silty clay loam 

0.089 

0.430 

0.010 

1.23 

1.68 

(10) Sandy clay 

0.100 

0.380 

0.027 

1.23 

2.88 

(11) Silty clay 

0.070 

0.360 

0.005 

1.09 

0.48 

(12) Clay 

0.068 

0.380 

0.008 

1.09 

4.80 
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2.4 Modified Picard (chord slope) approximation: 

The classic method of Celia et. al. (1990) is based on a backward difference 
approximation to dOldt in first order ODE system: 

M— + K/i = F (2.7) 


n n+I ft 


This non-linear system can be solved using a modified fixed-point (modified 
Picard) iteration: 

^/H-I.m+Uy.n+l.m^j^/n+1 ^2.10) 

(2.9) and (2.10) will be referred to as the Celia et al. solution to Richards 
equation. In the Celia et al. scheme, d9 Idh in the C matrix is evaluated analytically. 
Rathf elder and Abriola (1994) showed that the Celia et al. scheme is equivalent to a 
pressure-based backward Euler formulation, with the specific capacity dO Idh 
approximated using a chord-slope estimate: 

deX*' _ ^9 9"*^ -9" 

dh] ~AA A""' -A" ^ 

This approximation, proposed by Celia et al (1990), has been used in formulation 
of flow in present study, which has been discussed in chapter 3. 
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Chapter 3 

Problem formulation and simuxation 


3.1 General 

The whole study in this dissertation is aimed at removing the cumbersome 
numerical procedure and difficult experimental data, for the identification of VGM 
parameters for different soils. Presently, most of the available methods are based on 
head-measurement and those, which are based on outflow, too require pressure head at 
some location of soil column. 

As for as experimental convenience is concerned, head-measurement requires 
sophisticated transducers and other accessories to maintain the required precision ol 
experimental data. On the other hand, discharge measurement is comparatively simple 
process and data can be easily obtained. 

So, the objective is to calibrate the VGM parameters and develop a procedure foi 
which the only information to be supplied, is outflow of soil column. This will enable us 
to extract suitable values of VGM parameters more easily, to simulate the unsaturated 
flow mathematically. 

The work of this dissertation have been done in two steps: 

1 . Simulation of V GM outflow for possible range of a and n . 

2. Analysis of various data sets obtained from previous step. 



Above two steps have been resulted in the shape of equations from which the 
parameters can be calculated by extracting two factors, namely, Tso and 56o-4o-, from given 
outflow data. These factors will be discussed in detail in the section 3.4.1. 

Before starting the discussion on outflow data, it is important to discuss the 
various forms of Richard’s equation and their applicability. 


3.2 Suitability of various forms of Richard’s equation 

Three forms of Richard’s Equation as given in chapter 1, can be written as below: 


dt dz 

cm- 

dt 

dt dz 


Did) 


dz 


dK(0) 

dz 


dz 


Kih) 


m 

dz 


dKjh) 

dz 


Kid) 


dz 


dKjh) 

dz 


(Moisture-based ^form) (3.1) 


(Pressure-based, A-form) (3.2) 


(Mixed or Coupled form) (3.3) 


The terms on left side of above equations describe the effects of draining and 
filling pores, so statement in terms of the temporal change in moisture content is more 
appropriate than description via pressure. In other words, the term i^&d}i)dhldt is more 
appropriately written in its simpler form d&dt. 

On the right side of equations, it is noted that the spatial derivative of the 
hydraulic head is used to describe the driving force for fluid movement. This is the most 
direct mathematical statement of the fact that head differences do indeed supply the 
energy required to move fluid. Specification of hydraulic conductivity as a fimction of the 
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pressure head, K{h), is, however, not directly representative of underlying physics. It is 
water-filled pores which facilitate transmission of water through the porous medium. 
Although, water content and pressure are directly related but formally, the hydraulic 
conductivity should be expressed in terms of the moisture content, K{0). Besides, Kih) 
typically exhibits hysteresis, where as in K{6) this phenomenon is generally less 
pronounced [Nielsen (1986)]. 

Hence, the most direct mathematical expression of the physics of unsaturated 
flow, given the consideration above, is the mixed form presented in equation (3.3). 

The present study, as explained earlier, intend to deal with discharge outflow and 
hence the boundary condition, considered here, is of Neumann type (constant flux) 
Boundary Condition, i.e., 

= 0<z<L ■■■ ■■■■ 

= t>0 

q{Q,t) = q^, t>Q 


3.3 Formulation 

As discussed in chapter 1, finite difference scheme is used in present study. The 
standard finite difference method (FDM) for mixed form of Richard’s equation obtained 
by a backward Euler method for temporal discretisation is 

A/ 

(3.4) 


K. 


Mn 




n-¥l,m 








Az 
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where AT,"*';" indicate the inter-block hydraulic conductivities. The key to the 
method is expansion of {n a truncated Taylor series with respect to h, about the 

expansion point h "*''" , namely. 


e: 




dh 




(3.5) 


If all terms higher than linear are neglected in eq. (3.5) and it substituted into 
eq.(3.4), it results 
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rt+l,/n 




f \ / 
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n-fUm+l 7 «+l.m+I 
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(3.6) 


Substituting the increment of pressure head at subsequent iteration levels, i.e.. 


from eq. (3.6), we get 
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(3.7) 


where is defined as the residual associated with the modified Picard 

iteration [Celia et al, (1990)]. 
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Now these equations (3.7) can be written in matrix form as 
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where the coefficients O, P, Q and R are: 


Qn+Um 


{AzY 


p; 




Cl 




•+- 


K, 


n+i,m 

1 - 1/2 


+- 


a:, 


1 + 1/2 


e; 


/f+l,/n 


R 


fti-Um 


At 

(Az) 

_ -^1+1/2 

and 

{Azf 


___ ‘^/+l/2 ll 

/r+l,m 

(Az)^ ‘‘ 

^•+1 “ 


(AzY 


rn-hl,m \ ^/-l / 2 

/ / a -\2 \^i 


n+l,m 


/ f + I,JW 


i^Y 



T^n+i,m 

^/+l/2 


-K 

Az 


M/2 


_ Qft 


A/ 


(3.8) 


20 



3.3.1 Discretisation of boundary condition 


In order to apply the boundary conditions eq. (3.8) is modified at the top and 
bottom node. In present study, Neumann type of boundary condition is applied, where 
flux is known at top node. It is given by 


q = -\ 


dz 


(3.9) 


where q is the prescribed flux. 


3.3.2 Estimating the inter-block hydraulic conductivity 

The application of FDM provokes the problem of approximating the i ± ‘/z label in 
space, referred to as “weighting” necessary for the determination of interblock hydraulic 
conductivity values Different weighting formulas for estimating interblock 

quantities from the available grid point are proposed {Havekamp and Vauclin (1979). 
Arithmetic mean 




"/±l/2 


Geometric mean 

Harmonic mean 




sr^ifccT A* 
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Upstream method 
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Srivastava and Guzman (1995) has analyzed various formulas and recommended 
the use of geometric mean. In present analysis, the Richard’s equation is approximated 
through modified Picard’s scheme and solved with finite difference method. Geometric 
mean is used to carry out the inter-block hydraulic conductivity (i^r,+i /2 and Ks.\a). 


3.4 Numerical simulation 


for the numerical simulation all the quantities are considered in dimensionless 
form as follows: 


h' =hlL 
K’ =KIK, 
9-e^ 


( 9 '=- 


a =a-L 
q=q/K, 
t-K, 


0 ,- 0 , 




L-{0,-9,) 


. Total depth of soil column jL = 1 00cm 
, Constant Darcy’s flux at top ^top = 0.5 

In the present study, the range of a is taken firom 0.02 to 0.15 cm’* and that of n 
is taken from 1.4 to 3.0. 
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Chapter 4 

Results and discussion 


4.1 Outflow vs. Time: Effect of a and n 

Figure 3.1 and 3.2 shows pattern of outflow variation when n and a has been 
varied respectively, keeping all other factors constant. 

Now, by observing the trends of these graphs, it can be easily noticed that the 
variation in the outflow, can mainly be described by two factors, namely 

1 . The time when the change in outflow occurs. 

2. How rapid the outflow is increasing at saturation. 



Fig. 4. 1 ; Pattern of outflow variation with n. 
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Fig. 4.2: Pattern of outflow variation with a varying. 


Defining two quantities, Tso and 56o-40, as follows, can mathematically, represent 
the above two factors. 

Time as T^n : 

It is the time when outflow becomes 50% of the constant inflow at the top of soil 
column. As the graph shifts towards right, T 50 will increase. Indeed, the time factor can 
be analysed by any T, not necessarily T 50 . Still, T 50 is better choice and this can be 
explained by the fact that for the estimation of parameters experimental data has to be 
used and hence the input ( i.e., T ) may have significant errors as it approaches to 
saturation. Error will also be dominating at the points where the outflow just starts 
increasing. In other words, the curved part of outflow may result in much higher degree 
of errors. Therefore, the best option is to be at comparatively straighter middle part. Thus 
Tso is a better choice of input data. 
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5 Iope as Sm-m : 


It is clear from the above graphs that the slope of curve is not constant. In any 
outflow curve for some particular value of a and n first the slope is increasing from zero 
to some particular value S and then again decreasing to zero. This maximum value of 
slope, i.e., S is of interest as it is varying smoothly with a and n. The slope thus described 
can represented as the difference of time T for two particular value of outflow: 

‘S'60-40~T60 — 740 (4. 1 ) 

where Jeo and r 4 o are defined in same fashion as Tsq. 

Again lot of choices of T, are there. As it has been explained earlier that curved 
part should be avoided, but here the interesting aspect is that the difference must be large 
enough to overshadow the possible error in the experimental data. After analyzing 
various combinations slope defined in eq. (4.1) has been found to be the better choice. 

FDM equations discussed in Section 3.3, have been simulated in dimensionless 
form, therefore, making Z 50 and 56 o- 4 o a dimensionless quantity. In the next step of 
present work, these two quantities have been studied and used for determining the VGM 
parameters. 

Since the objective of this study is to estimate the parameters from outflow, only 
the outflow has been studied, but similar trends can be noticed in the pressure head. 
Figure 3.3 and 3.4 shows the pressure head distribution at steady state for changing 
values of nr and n. 
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4.2 Variation of Tgo and S^q^q 

Fig (4.5) and (4.6) shows the pattern of T 50 with respect to a and n respectively. It 
is very clear from the graphs that T 50 is increasing with both a and n. the sensitivity of Tso 
is more w. r. t. a but in fig 4.5, as n increases the sensitivity is decreasing. 

The curves obtained in fig 4.5, has been fitted in mathematical expression given 
in eq. (4.2) in order to obtain T^o, if the two parameters are known. 

(4.2) 

where coefficients a and b, both are dependent of n. 



0.02 0.04 0.06 0.08 0.1 0.12 0.14 

a incnT^ 

Fig. 4.5: variation of T 50 with a 
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Fig. 4.6: variation of T 50 with n 


Various values of coefficients of eq. (4.2) are given below in Table 4.1 


Table 4.1 Values of coefficients of eq. (4.2) for different n 


a 

B 

1.05069 

17.93729 

1.18563 

18.40056 

1.2902 

18.86338 


1.37096 


1.43307 


1.48047 


1.51634 


1.54333 


1.56343 


19.34075 


19.84663 


20.38249 


20.9446 


21.52271 


22.11226 






























2.3 

1.57817 

22.70734 

2.4 

1.58877 

23.3024 

2.5 

1.59616 

23.89382 

2.6 

1.60103 

24.47887 

2.7 

1.60395 

25.05784 

2.8 

1.60533 

25.6258 

2.9 

1.60557 

26.18049 

3.0 

1.60489 

26.72085 


Again, these coefficients are fitted in polynomial models given by following 
equations 


a 


1.2355 




n 


6.2945 


NO.. 


+ 1.61946 


1 + 

^ 1.36362 


(4.3) 


and 


6 = 9.85517 + 5.6091 •« 


(4.4) 


Now, eq. (4.2)-(4.4) can give the value of Tso for any value of O' and n 
Similarly, the expression is also obtained for Sso+o- Fig (^-7) ^nd (4.8) shows the 
pattern of with respect to a and n respectively. Here, the point to be noticed is that 
the slope increases with increasing n, but decreasing with an increase in a. 
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The various curves for 560-40, obtained in fig 4.7, has been fitted in mathematical 
expression given in eq. (4.5) in a similar fashion as in T50, if the two parameters are 
known. 

^60 ~T^q=T^-{- A- e~°'‘ (4.5) 

where 

Tc =—0.0117 + 0.03682 -/i- 0.01 136 -Ai* +0.00122-/j^ (4.6) 

= -0.334 + 0.41776-/j-0.1459-/i- +0.018-/i^ (4.7) 

and, 

t = 0.0531 + 0.09482- « + 0.0337 •«- + 0.0039 -ai^ (4.8) 

Table 4.2 shows various values of these coefficients for various n. 


Table 4.2: variation of coefficients of eq. (4.5) w.r.t. n 


n 

Tc 

A 

t 


0.00158 

0.0136 

0.02482 


0.00258 

0.02543 

0.02572 

1.6 

0.00353 

0.03499 

0.02844 

1.7 

0.00451 

0.04338 

0.03 

1.8 

0.00549 

0.05061 

0.03107 

1.9 

0.00626 

0.05643 

0.03229 

2 

0.00687 

0.06151 

0.03317 

2.1 

0.00741 

0.06603 

0.03369 

2.2 

0.00786 

0.06998 

0.03405 

2.3 

0.00828 

0.07346 

0.03426 

2.4 

0.00868 

0.07673 

0.03422 

2.5 

0.00901 

0.07968 

0.03415 
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2.6 

0.00926 

0.0823 

0.03419 

2.7 

0.00946 

0.08465 

0.03422 

2.8 

0.00972 

0.08698 

0.03397 

2.9 

1 0.00993 

0.08914 

0.03371 

3 

0.01011 

0.09112 

0.03347 


4.3 Relationship of Tso and S^o^o with parameters 

Relationship of parameters and outflow data can be presented in two forms, 
graphical and mathematical. But the problem in determining mathematical form persist in 
the fact that this relationship will be based on the regression analysis of data obtained 
from the fitted equations for T 50 and ^’ 50-10 [eq. (4.2) and (4.5)]. This results in increased 
errors of final output, i.e., estimate of parameters. 


4.3.1 mathematical expression for parameters 

After having equations (4.2) and (4.5) values of a and S 60-40 can be obtained while 
changing Tso and n as independent variables. Once again repeating the similar process of 
curve fitting, expressions for a and n have been obtained as follows. 


^ — U + 6 • iS'50_4^ C • + d ' >5'6o-40 

(4.9) 

where 


fl = 1.34853- 0.47317 -r^o + 0.74575 • -0.32916-r/o 

(4.10) 

b = 18.5956 + 1.12316 ■ 

(4.11) 


(4.12) 

r/ = 1.03x10" 

(4.13) 
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Once n value is obtained, a can be estimated from eq. (4.2) as. 


b \ a) 

where a and b can be obtained from eq. (4.3) and (4.4) respectively. 


( 4 . 14 ) 



Fig. 4.9: Error in a and n back-calculated from eq (4.9) & (4. 14), for 252 data sets 

obtained from simulation 

Fig. 4.9 shows the % error in parameter estimate from eq. (4.9) and (4.14). As the 
errors are going beyond acceptable limits, the mathematical expression derived above, 
are only reliable for a very short range, a =0.03 to 0.14 and « = 1.9 to 3. it is clear from 
the range that equations are valid only for coarser soils. 
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The main problem persist with the repeated regression, and hence the value of 
parameters is highly sensitive to the coefficients of eq. (4.9) and (4.14). Even a very 
small variation in these coefficients can severely affect the output. 

4 . 3.2 Graphical representation of Tso and .^60-40 to estimate parameters 

In view of above discussion, it is clear that a well-correlated relationship can be 
established between outflow and flow parameters. Fig 4.9 illustrates a graphical 
representation of Tso and .$' 60-40 showing the influence of parameters. 



Fig. 4.10; 560-10 vs. fso : graphical representation of studied domain of numerical values of 

parameters 
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Fig. 4.11: Tso and 5'6o.40 for higher values of a and n 

This chart showing various combinations of T 50 and <S' 6 o -40 can easily be used for 
estimation of VGM parameters. At the higher values of a and « fig. 4.8 shows very 
congested points, hence this part is shown in fig 4.10 for convenience. 

Now, using fig. (4.9) and (4.10) values of parameters can be estimated, knowing 

the Tso and 560-40 from outflow data. 
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Chapter 5 

Conclusions and future study 


5.1 Conclusions 

The objective of this dissertation is to study the outflow in order to give a suitable 
method for estimation of flow parameters. As discussed in Chapter 4, the graphical 
representation of outflow quantities {Tsq and -Seo-io) is quite convenient way to get an 
estimate of or and n. The whole study can be concluded in following points: 

• Outflow can be very smoothly correlated with the unsaturated flow parameters, by 
defining the two quantities, Tso and 560-40- 

• Any mathematical model is difficult to obtain for the whole range of a and n due 
to two reasons: 

(1) Both parameters are highly sensitive to the input. Even for a small 
error in input (Tso and 56o-4o) can lead to an unacceptable error in 
parameter estimate. 

(2) Since the input has to be an experimental data, it caimot be expected to 
have perfect precision. 

• A graphical representation has been given, from which the parameters can be 
easily estimated, knowing Tso and 56o-4o- 
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5.2 Recommendations for future work 


Present study is mainly focused on outflow data and the emphasis is on soils of 
sandy nature (the smaller values of a and n have not been analyzed). The phenomenon 
may be quite different for very small values of a and n. For smaller values of parameters 
(a < 0.02 and n < 1.4) the convergence is very difficult to obtain. So, with the 
development of a better convergence scheme, the phenomenon should be studied for very 
small values of parameters. 

The boundary condition, itself, can make significant changes in the outflow and 
thus in Fso and Hence the effect of boundary condition needs to be studied. 

In the proposed method main focus is on outflow data, but the pressure head (as 
discussed in Chapter-4) can also be used to estimate the parameters. 
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